package scu.maqiang.elasticity;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES3H82;
import scu.maqiang.mesh.Mesh3H8;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.NewIterSSolver;
import scu.maqiang.numeric.SRMatrix;

public class ElasticPressureLoad3DH82 {

	public static void main(String[] args) {
		long startTime = System.nanoTime();
	    double E = 2.0e11, Nu = 0.3e0;
		double verbosity = 1;

		Mesh3H8 mesh3D = new Mesh3H8().cube(80, 10, 10);
		mesh3D.scale(8.0e1, 1.0e1, 1.0e1);
		ScalarFunc dispRegion = (xyz, label, param) -> Math.abs(xyz[0]) < 0.00001 ? 1.0:-1.0;
		int dispLabel = 10;
		mesh3D.setBoundaryLabel(dispRegion, null, 10);
		
		ScalarFunc forceRegion = (xyz, label, param) -> Math.abs(xyz[0]) > 79.99999? 1.0:-1.0;
		int forceLabel = 20;
		mesh3D.setBoundaryLabel(forceRegion, null, forceLabel);
		
		FES3H82 fs= new FES3H82(mesh3D);
		
		SRMatrix K = new SRMatrix(fs.GetNdof());
		double[] x = new double[fs.GetNdof()];
		double[] RHS = new double[fs.GetNdof()];
		double[][] disp = new double[3][fs.GetNdof() / 3];
		double[][] strain = new double[6][fs.GetNdof() / 3];
		double[][] stress = new double[6][fs.GetNdof() / 3];
		
		double[] elasticCoef = {E, Nu};
		fs.assembleStiff(elasticCoef, BVPType.COMMON, K);
//		K.sort();
//		System.out.println(K);
//		System.exit(0);
		double[] surForce = {0.0, 0.0, -1000.0};
		fs.assembleSource(surForce, BVPType.COMMON, RHS);
		//fs.assembleSurfaceForce(surForce, BVPType.COMMON, RHS, forceLabel);
//		System.out.println(MVO.toString(RHS));
//		System.exit(0);
		//施加边界条件
		//....
//	    label(1) = dispLabel
//	    !call fs.applyBCZero(K, RHS, label)
		fs.applyBC_MBN(K, RHS, Direct.All, 0.0, dispLabel);
//	    call fs.applyBC(K, RHS, label, 0.0D0, 0.0D0, 0.0D0, 4)
		NewIterSSolver solver = new NewIterSSolver(K);
//	    call solver.initIterSSolver
//	    !call solver.PCGSSOR(K, RHS, x, 1.5D0, 2)
//	    call MKL_DSS_Solve(K, RHS, x)
		solver.PCGSSOR(RHS, x, 1.5, 1);
		
		fs.extractComponent(x, disp);
	    System.out.println("Displacement:");
	    System.out.println("u min : " + MVO.min(disp[0]) + " max: " + MVO.max(disp[0]));
	    System.out.println("u min : " + MVO.min(disp[1]) + " max: " + MVO.max(disp[1]));
	    System.out.println("u min : " + MVO.min(disp[2]) + " max: " + MVO.max(disp[2]));
		
	    mesh3D.toTecplot("cantilevel3DH82Sol_1.dat", disp);
	    mesh3D.toTecplotDeform("deformCantilevel.dat", disp, 10.0);
//	    mesh3D.toTecplot("cantilevel3DH82Sol_1.dat", strain);
	    
	    
//	    label(1) = dispLabel
//	    !call fs.applyBCZero(K, RHS, label)
//	    call fs.applyBC(K, RHS, label, 0.0D0, 0.0D0, 0.0D0, 4)
//	    call solver.initIterSSolver
//	    !call solver.PCGSSOR(K, RHS, x, 1.5D0, 2)
//	    call MKL_DSS_Solve(K, RHS, x)
//	    
//	    uvw = sqrt(u**2 + v**2 + w**2)
//	    write(*, *) "Displacement Vector Sum:"
//	    write(*, *) "min : ", minval(uvw), "max: ", maxval(uvw)
//
//	    allocate(ex(fs%nDof/3), ey(fs%nDof/3), ez(fs%nDof/3))
//	    allocate(exy(fs%nDof/3), eyz(fs%nDof/3), exz(fs%nDof/3))
//
//	    call fs.computeStrain(u, v, w, ex, ey, ez, exy, eyz, exz)
//	    write(*, *) "Strain:"
//	    write(*, *) "ex min : ", minval(ex), "max: ", maxval(ex)
//	    write(*, *) "ey min : ", minval(ey), "max: ", maxval(ey)
//	    write(*, *) "ez min : ", minval(ez), "max: ", maxval(ez)
//	    write(*, *) "exy min : ", minval(exy), "max: ", maxval(exy)
//	    write(*, *) "eyz min : ", minval(eyz), "max: ", maxval(eyz)
//	    write(*, *) "exz min : ", minval(exz), "max: ", maxval(exz)
//
//	    allocate(tx(fs%nDof/3), ty(fs%nDof/3), tz(fs%nDof/3))
//	    allocate(txy(fs%nDof/3), tyz(fs%nDof/3), txz(fs%nDof/3))
//
//	    call fs.computeStress(ex, ey, ez, exy, eyz, exz, E, Nu, tx, ty, tz, txy, tyz, txz)
//	    write(*, *) "Stress:"
//	    write(*, *) "tx min : ", minval(tx), "max: ", maxval(tx)
//	    write(*, *) "ty min : ", minval(ty), "max: ", maxval(ty)
//	    write(*, *) "tz min : ", minval(tz), "max: ", maxval(tz)
//	    write(*, *) "txy min : ", minval(txy), "max: ", maxval(txy)
//	    write(*, *) "tyz min : ", minval(tyz), "max: ", maxval(tyz)
//	    write(*, *) "txz min : ", minval(txz), "max: ", maxval(txz)
//
//	    allocate(eqStress(fs%nDof/3))
//	    call fs.computeVonMisesStress(tx, ty, tz, txy, tyz, txz, eqStress)
//	    write(*, *) "Von Mises Stress"
//	    write(*, *) "min : ", minval(eqStress), "max: ", maxval(eqStress)
//
//	    call mesh3D.exportSolution("cantilevel3DH82Sol_1.dat", u, v, w)
//	    call mesh3D.exportSolution("cantilevel3DH82Sol_2.dat", ex, ey, ez, exy, eyz, exz)
//	    call mesh3D.exportMeshDeform("cantilevel3DH8Deform.dat", u, v, w, 5.0D0)
	    
        long endTime = System.nanoTime();
        System.out.println("Elapse Time: " + (endTime - startTime)/1e9 + " s");
	}

}
